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Abstract. We obtain a fourth order accurate numerical algorithm to integrate the 
Zerilli and Regge-Wheeler wave equations, describing perturbations of nonrotating 
black holes, with source terms due to an orbiting particle. Those source terms contain 
the Dirac’s delta and its first derivative. We also re-derive the source of the Zerilli and 
Regge-Wheeler equations for more convenient definitions of the waveforms, that allow 
direct metric reconstruction (in the Regge-Wheeler gauge). 
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1. Introduction 

In the special case of perturbations around an spherically symmetric background, as 
the Schwarzschild one with mass M, we can benefit from the decomposition into 


spherical harmonics (labeled by m) of the metric and hence obtain a set of Hilbert- 
Einstein HU EH held equations for metric coefficients depending only on time and 


radial coordinates, i.e. the Schwarzschild’s ( t,r ). Zerilli has written down the 
General Relativity held equations in this decomposition, for the Regge-Wheeler gauge 
[See Ref. 121| where some misprints to this equations have been corrected.] From these 
equations one can deduce wave equations for the two polarizations of the gravitational 
held, denoted by the waveforms odd ) ■ 

Our task is next to numerically integrate the Zerilli jSEj and Regge-Wheeler EH 
equations for even and odd parity perturbations respectively 


dl - a 2 - v, (even ’ odd) 



(i) 


= m^en,o dd) d r S[r - R(t )] + G(t)g en , odd} 5[r - R(t)}, 


where r* = r + 2Mln(r/2M — 1) is a characteristic coordinate, and 

[A 2 (A + l)r 3 + 3A 2 Mr 2 + 9A M 2 r + 9 M 3 ] 



( 2 ) 
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v (° dd ) _ 2 


1 - 


2 M 


A + 1 


-3- 


M 


( 3 ) 


are the Zerilli’s and Regge-Wheeler potentials respectively, and where A = (i — l)(l + 

2 )/ 2 . 

We have previously studied the lieadon collision of binary black holes in the extreme 
mass ratio regime|22, 230, E] and have developed and algorithm to deal in the time 
domain with the Dirac-delta (and derivative of) source term appearing in the Zerilli wave 
equation. This technique applies to both, even and odd parity perturbation equations 
of a Schwarzschild black hole. Integration of these perturbation equations j2>] in the 
time domain is, in general, much more efficient than in the frequency domain |22j . and 
can be directly related to full numerical simulations U HO El HUB 

There are at least two main motivations to seek for a higher than second order 
convergent numerical algorithm. The first comes from the need to compute radiation 
reaction corrections to the orbital motion of a particle. In order to do that we need to 
compute the ’self force’ that involves third order derivatives of the waveforms (in the 
Regge-Wheeler gauge jlSj) at the location of the point-like particle pEJ. The second 
important motivation comes from the need to have an efficient algorithm to generate the 
huge bank of templates needed for the analysis of data coming from the gravitational 
wave detectors such as LIGO and LISA ESI- Fourth order full numerical relativity 
has recently been implemented JEU] bearing in mind the same motivations. Also a 
fourth order numerical algorithm has recently been developed to deal with vacuum 
perturbations of Kerr black holes j3Ej- For the sake of completeness let us note that 
incorporating mesh refinement techniques can further help on the aspect of speeding up 
the generation of templates. 

The first order metric perturbations techniques of an spherically symmetric system 
split the fields into the two decoupled polarizations: For even parity perturbations we 
consider the following waveform 123] in terms of generic metric perturbations in the 
Regge-Wheeler notation [35] 

r — 2M ( H tm ^ rdrK tmj 


«n(r,f) = 


+ 


(A + 1) L 

r — 2 M 


K lm + 


A r + 3 M 
[r 2 d r G lm - 2h[ m \ , 


A r + 3 M 

This is related to Zerilli’s [33 even parity waveforms f>z er e 


by 




im 

Zer,even 


= d t ip, 


im 


+ 


47 t i \[21 


2 M) A 


(i) 

im 


(A + 1) (Ar + 3 M) 


where for an orbiting particle 


411 = irnoVz 6 l r ~ 


and it relates to Moncricf’s [ 22 ] waveform 'fff on even by a normalization factor 


ib im = 

t even 




im 

Mon,even 


(A + 1) 


( 4 ) 


( 5 ) 


( 6 ) 


( 7 ) 
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The contribution of the even modes represented by our waveform to the total 
radiated energy (either to infinity or onto the horizon) is given by 


— = — V ^ + 2 ) ! (d il> lm ) 
dt 64vr^ U~2)\ V Weven) 

im v ' 


( 8 ) 


For instance, for circular-equatorial orbits the source term of the even parity wave 
equation (JTJ) takes the simple form 


F'evenif) 


m 0 U°{R-2M) 3 - 

(A + 1) (A R + 3 M) R 2 Ylm 


(9) 


Gevenif) 8 


nm 0 (R-2M) U° (|<F) 2 (m 2 - A - 1 )Y tm 
A (A + 1) 


Q Y £m 7T7n 0 U° (i?-2M) 2 (!$) 2 
+ (A + 1) (XR + 3 M) 

0 Y irn nm 0 U°(R-2 M f (R 2 A + R 2 A 2 + 15 M 2 + 6 A RM) 

— 8--, (10) 

R 3 (A + 1) (XR + 3 M) 2 

where Yi m are the usual spherical harmonics, an overbar represents complex conjugation, 
and R, 0, $ define the trajectory of the orbiting particle in spherical coordinates. For 
the general orbit form of the source terms see | Appendix A[ 

The Regge-Wheeler wave equation describes the odd parity modes. We will consider 
the following waveform in terms of generic metric perturbations in the Regge-Wheeler 
notation 


«M) = j 


r 2 d r 


h £m [ 


r, t) 


dth l ™{r , t) 


Or _ 

= -v / l-2M/r/^”.(ll) 


This waveform is related to the Zerilli’sPHIJ and Moncrief \s |29| odd parity waveforms 




im _ nUim 


Zer,odd 


= if 


Mon,odd 


,lm _ (! - 2M / r ) 

YZer.odd. ~ 

ry 



by (See Eq. (1A.15I1 1 


d t if lm 


2 YTer^dd 


87iir(r -2M)Q em 

XVX+1 


( 12 ) 


(13) 


and to the Cunningham et al. m waveform ipQ 1 by 


if £m 


-2 


(1-2)! 

(e + 2 )! 


Yh m = 


1 V'g" 

2 A(A + 1) ’ 


(14) 


And very close to the Weyl scalar Im iff" 1 = [Here we used the 

Kinnersley tetrad, in the Schwarzschild background, and decomposed Im'k 2 into 
spherical harmonics]. 
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One of the advantage of this odd parity waveform definition over Zerilli’s is that it 
allows an straightforward construction of metric coefficients in terms of the waveform 
(and its time derivatives) in the time domain (see Ref. [23]). It is also possible 
to construct metric coefficients from the Zerilli waveforms, but it involves extrinsic 
curvature perturbations J3]. 

The contribution of the odd modes to the total radiated energy is 


dE 

dt 


1 (f + 2)! 
64 ^^ (£ — 2 )! 

im v ' 




(15) 


For circular, equatorial orbits the source term of the odd parity wave equation 
takes the form 


F(t) = 8 Tm °f ,? )(B : 2M>2 ( dgY em (Q, #), 


G(t) = 


A(A + 1 )R 


8wm 0 U°(t) (R - 2M) /dd>\ ^ 


\ dt J 

fdslE 

\dt J 


(16) 

(17) 


A(A + l)fl 

and for the general form of the source terms see Appendix A| 

For the sake of notational simplicity, from now on we will drop the (£ m) multipole 
and ’even/odd’ indexing from and the potential V. 


2. Numerical Method 


In this section we describe the algorithm used to integrate the wave equation CD 
numerically. While the left hand side of this equation is straightforward to integrate, 
the source S im contains terms with the Dirac’s delta and its derivative. Since we have 
not found 23] hi the literature a discussion of the numerical treatment of such sources, 
we shall describe the method in some detail. 

It is convenient to use a numerical scheme with step sizes A t = | Ar* = h, and with 
a staggered grid. This represents a characteristic grid. As Fig. Q] shows, this method 
connects points along lines of constant “retarded time” u = t — r* and “advanced time” 
v = t + r*. On this grid we have implemented a finite difference algorithm for evolving if 
with essentially no errors for integrating the wave operator beyond those due to round 
off. 

To derive our difference scheme we start by integrating Eq. Q over a cell of our 
numerical grid. Shown in Fig. [T] is the cell crossed by the orbiting particle (one of the 
four entering/exiting possible cases). We use the notation 



Applied to the derivative terms in Eq. © this gives: 



cL4 {-d t d t ip + = dA {-Ad u d v ip} = 









Fourth order accurate integration of black hole perturbations 


5 



tO+h 


to 


tO-h 


xO-h xO xO+h 


Figure 1. The cell that the particle crosses and the four areas defined to construct 
the second order algorithm. 


— 4 [fj(t + h , r*) + — h, r*) — if(t, r* + h) — ijj(t, r* — h)} . (19) 


Note that this result is exact; it contains no truncation errors. 

For cells through which the worldline passes, the integral of the source term in Eq. 
m must be evaluated. The integration over the cell, when done with due regard to the 
boundary terms generated by the 5'[r — R(t)], gives 


SdA 



G(t,R(t )) 

1 - 2 M/R(t) 


( F(t,r) \ 


VI - 2 M/r) 

r=R(t) 


±2 Fjt^Rjh)) 

(1 - 2M/R{t l )f 
± 2 F(t 2 ,R(t 2 )) 

(1 - 2 M/R(t 2 )f 


(l=Ffl*(«i)) 

(l ±&(h)) 


( 20 ) 


The f dt integral in the first term can be performed to any desired precision since F and 
G are known (analytic) functions. This integration can be considered as exact as well, 
since we assume the trajectory of the particle is known a priori , so we can perform the 
integral over the trajectory as accurate as needed by evaluating the F and G on as many 
points as we want (do not need to be grid-points). For our goal of quartic convergence, 
a Simpson approximation j33J for the integration is adequate. In the second term the 
upper (lower) sign is for particles entering the cell from the right (left), or equivalently 
for r* > r* (rjj 1 < r*). In the same way, in the third term the upper (lower) sign is for 
particles leaving the cell to the right (left), or equivalently r% > r* (rf < r*). 


2.1. Second order method 


We next consider the integration of the potential term over the cell, 
not crossed by the particle, then we can use a trapezoidal rule 



dA{-Vf>} 


-h 2 [V(r*)if(t + h, r*) + V{r*)i)(t - 


If the cell is one 


h, r*)+ 
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V(r* + h)if(t, r* + fi) + V(r* - h)if(t, r* - fi)] + 0 (h 4 ) . (21) 

The h 4 order error in a generic cell is equivalent to an overall 0(h 2 ) error in if after a 
finite time of evolution. 

The result in Eq. m assumes that if is smooth in the grid cell. It cannot be applied 
to those cells through which the particle worldline passes, since if is discontinuous across 
the worldline. For such cells we first obtain the coordinates (r*,£i) of the point where 
the particle enters the cell and (r^, t 2 ) where the particle leaves it (see Fig. |T| for one 
of the four entering/exiting possible cases). Next, we divide the total area of the cell, 
(4 h 2 ), into four subareas defined as follows: A 2 is the part of the area of the diamond 
below t — fi, A 3 is the part of the area of the diamond over t = t 2 , A 4 is the remaining 
area to the left of the particle’s trajectory, and A 4 is the remaining area to the right. 

The integral of the Vif term over the area of the cell is approximated by the sum 
of this function evaluated on the corners of the cell multiplied by the corresponding 
sub-area A*. This gives us 



dA{—Vif} = — V(r*) [if(t + fi,r*)A 3 + if(t,r* + fi)A 4 

+ if(t , r* - h)A 1 + if(t - h , r*)A 2 \ + 0(h 3 ). (22) 

The truncation error in each such cell is of order (fi 3 ), just enough to have quadratic 
convergence, since only one cell with the particle has to be evaluated per time step. 

I 11 summary, our numerical scheme, for cells through which the particle worldline 
does not pass, is to solve for if{t + h,r*), using Eq. (fill and Eq. m representing the 
integral over a single cell of the homogeneous version of Eq. ©• For cells containing the 
worldline, Eq. (fTH . Eq. dH, and Eq. (HU) are used instead. The evolution algorithm 
we use is then 


if(t + h,r*) 


—if(t — fi, r*) + [if(t, r* + fi) + if(t, r* 


fi)] 



+ £>(fi 4 ) , (23) 


for cells not crossed by the particle, and 

if(t + fi, r*) = — if(t — fi, r*) 

+ if(t,r* + fi) 

+ if(t, r* — fi) 


V (r*) 

1 + ^ 1 (a 2 -a 3 ) 


1 - 

1 - 


4 

V(r*) 

4 

V(r*) 


V (V*) 



[A 4 + A3) 

{Ai + A 3 ) 

sdA + e>(fi 3 ), 


(24) 


for the cells that the particle does cross [231123 

The above equations cannot, however, be used to initiate the evolution off the first 
hypersurface. If t = 0 denotes the time at which we have the initial data in the form of 
if(Q,r*) and d t if( 0 ,r*), we lack the values if(t = —fi), necessary to apply the evolution 
algorithm. We can, however, use a Taylor expansion to write 


if(—h , r*) = if(+h, r*) — 2hd t if(0, r*) + 0(h 3 ) . 


( 25 ) 
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The right hand side can be used in place of r*) in the the algorithm to evolve off 

the first hypersurface and algebraically solve for ^(+/q r*). It is important to note that 
this substitution is valid only if if(t,r*) is not singular between t = —h and t = +h. 
This requires that the particle worldline does not cross the vertical line at r* between 
t = —h and t = +h. In setting up the computational grid, we have been careful always 
to avoid such a crossing. 

For the special case when the source term does not contain derivatives of the Dirac’s 
Delta, i.e. F = 0, the waveform is continuous, i.e. C°. In this case we can still use 
the expression m for the integration over the potential term for all cells and obtain 
an overall error 0(h 3 ), thus producing a cubic convergent algorithm. This has a direct 
application for the integration of the Hilbert-Einstein field equations in the perturbative 
regime: In the harmonic gauge the general relativistic equations take the form of wave 
operators with source terms, given by the energy-stress-tensor, T^, proportional to 
Dirac’s deltas only (no derivatives of it) [Barack and Lousto in preparation.] 

3. 4th order Algorithm 

As we have seen the only part of the integration algorithm that needs to be approximated 
on the numerical grid is that of the potential term V(r*) i/j(t,r*). There are two cases 
to be considered separately: Those cells that the particle crosses and those that do not. 
For the later the integration is simpler and we will dealt with it first. 

3.1. Vacuum case 

For the sake of notational simplicity let us denote the integrand by 

g(t, r*) = V £ even ’ odd (r*) f)™™' odd (t, r*). (26) 

The integration over the cell of 



can be performed by a double Simpson [HI] integral providing the required forth order 
accurate evolutions. In order to perform this integral with the Simpson method we need 
to evaluate g at all the points shown in Fig. [3 The final result is 

2 

{#1 + 92 + #3 + #4 + 4 (#12 + $24 + #34 + #13) + 16#o} + O ) • 

(28) 

The use of this algorithm as it stands is possible but it implies to double the grid 
points in the time and space directions, and to retain information about four time slices, 
thus quadruplicating the number of points needed to evolve the same time interval. 
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xO-3h x0-2h xO-h xO xO+h x()+2h x()+3h 


Figure 2. Vacuum case: the cells that the particle do no cross. 


It is possible though to evaluate the extra points (marked as crosses in Fig. [21) from 
the original grid points (marked as circles in Fig.H. First, in order to avoid the central 
grid point, we can evaluate g 0 on the slice as follows 

go = 77 [9V^(x 0 - h) if (to, x 0 - h) + W^xq + h) if(t Q , x 0 + h) 

lo 

- Ve(x Q - 3 h) if (t 0 , x 0 - 3 h) - V e (x 0 + 3 h) if(t 0 , x 0 + 3 h)] 

+ O (h 4 ) (centered), (29) 

0 o = 77 [5V e (x A ) if (t 0 , x A ) + 15Ve(x A + 2 h) if(t 0 , x A + 2 h) 

lo 

- hV e (x A + 4 h) if (to, x A + Ah) + V e (x A + 6 h) if (to, x A + 6 h)] 

+ O (hf) (x A : left boundary), (30) 

9 o = TX [5 Vi(x B ) if (to, X B ) + 15Vg(x s - 2 h) if(t 0 , x B - 2 h) 

lo 

- 5 V e (x B ~ 4 h) if (t 0 , x B - Ah) + V e (x B - 6 h) if (to, x B - 6 h)} 

+ O (/r 3 ) (x B : right boundary), (31) 

Since the boundary values are only used once per time step, O (h 3 ) is all we need to 
evaluate g 0 . 

I 11 order to compute the points between the time levels of the original grid we use 
the second order evolution scheme, Eq. (El adapted to the current case 


0i3 + 012 = V e (x 0 - h/2) (if 13 + ifi 2) 

= V e (x 0 - h/ 2 ) (ifi +if 0 ) 

024 + 034 = Vz(Xo + h/2) (if 24 + ^ 34 ) 

= Ve(x 0 + h/ 2 ) (ifo + ifi) 1 - 


1 fh 


2 V 2 


1 fh 


2 \ 2 


V/>(x 0 - h/ 2 ) 


Ve(x 0 + h/2) 


(32) 
+ O {h 4 ) , 

(33) 

+ O ( h 4 ) . 


This completes the computation for the cells not crossed by the particle. Note that 
including two more points (the circles labeled as Xo ± 3 h in Fig. H) in the algorithm 
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allowed us to avoid the evaluation of the held in five new points (the crosses in Fig. |2J)). 


3.2. Cells with the particle 


The key idea here is to expand the function g(t,x = r*) = V(r*) ip(t,r*) around the 
left/right vertexes of the cell labeled with 1 and 4 in Fig. [21 and integrate the left and 
right parts of the cell as defined by the trajectory of the particle. We next write down 
explicitly the four possible integrals to the left and the four to the right depending on the 
trajectory of the particle (entering/exiting the cell). Here we will assume the function 
g(t, x ) can be expanded as individual Taylor series on both sides of the of the trajectory 
of the particle. So, to the required order (since we only use this algorithm once per time 
step) we can write 


dg 

g R , L (t, x) = g[t 0 , x 0 ±h} + ^[t 0 , x 0 ±h] (x-x 0 t h) + 

dg r , , , i ,. . , . d 2 g (x-x 0 Th) 2 

— [to, X 0 ± h\ (t - to) + ^2 [to, X 0 ± h\ - - -h 


d 2 g 

dtdx 
}2 


[to, xq ± h] (t - to) (x - x 0 =f h) + 


^ Kxo±h] (lxJsL + om . 


(34) 


We give the explicit construction of the derivatives of the function g(t, x) out of the 
evaluations of g(t,x ) at nearby grid points in the Appendix B 


3.2.1. Left side integral : 

i) This case is displayed in Fig. Eland the integral over the potential term (in the 
shaded area) reads 


rt 2 rx P (t) 


gL^t , x) dt dx 


J XQ —h 
rto 

dt 


rx p (t) 


gL(t,x ) dx + 


r*2 


rx p (t) 


dt 


gi,(t, x) dx. 


' VQ—t — h 


'to 


' t — h — UQ 


(35) 


Here t\ and t 2 are respectively the times the particle enter and leaves the cell. Its radial 
trajectory given by x p (t). Note that the central point of the diamond (u 0 = t 0 — x 0 , v 0 = 
t o + Xo) is not a grid point. 

ii) This case is displayed in Fig. E] and the integral over the potential (in the shaded 
area) term reads 


Hi r x p(t) 
J tQ — h J xq — h 
+ 


f 

'to—h 


r*t-\-h—uo 


g R (t,x) dtdx— j dt gi,(t,x)dx 

VQ—t—h 

Ho r x pd ) H 2 rx P {t) 

i dt / gi.(t, x) dx + / dt gi(t,x)dx. 

t\ J VQ — t—h J tQ J t—h—UQ 


( 36 ) 








xO-h xO xO+h 

Figure 3. Case i), the particle enters (at time t \) and leaves (at time t?) the cell on 
the left side. 


xn(t) 



tO+h 


to 


tO-h 


xO-h xO xO+h 


Figure 4. Case ii), the particle enters (at time ti) on the right and leaves (at time 
£ 2 ) the cell on the left side. 


iii) This case is displayed in Fig. 0 and the integral over the potential term (in the 
shaded area) reads 


cto+h rx p (t) 


gL(t, x ) dt dx = 


' xq— h 


rt 2 


r x p (t) 


+ 


dt 


'to 


rto rx P (t) 

dt / gi(t,x)dx 

J VQ—t—h 

pto+h nvo —t+h 

I gi(t,x)dx. 

t — h — UQ 


dt 


(37) 


'*2 


/ g L (t,x)dx + 

' t—h — UQ 

iv) This case is displayed in Fig. IHland the integral over the potential term (in the 
shaded area) reads 


rto+h rx p (t) 


g L (t^x) dtdx = 


rti 


rt+h—uo 


dt 


'to—h J xq—H 


'to-h 


VQ—t—h 


g L (t,x) dx 
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Figure 5. Case iii), the particle enters (at time t\) on the left and leaves (at time tf) 
the cell on the right side. 



xO-h xO xO+h 


Figure 6. Case iv), the particle enters (at time ti) and leaves (at time £ 2 ) the cell on 
the right side. 


+ 

+ 


rto rx P (t) 

/ dt g L (t, x) dx + 

' t\ JVQ—t—h J 

rto+h f'VQ—t+h 

/ dt g L (t,x)dx. 

'$2 J t — h — UQ 


rt2 


to 


dt 


rx p (t) 


' t—h—uo 


g L (t,x)dx 


(38) 


3.2.2. Right side integral : 

i) This case is displayed in Fig. Eland the integral over the potential term (in the 
11011 -shaded area) reads 


rto+h rXQ+h. 


'to—h J x p (t) 


g R (t, x) dt dx = 


dt 


rt 0 


pt-\-h — UQ 


+ 


dt 




' Xp(i) 


gn(t, x) dx + 


’to~h 
rt 2 


rt-\-h -uq 

' VQ—t—h 
pvq— t-\-h 


g R (t,x) dx 


dt 


'to 


Xp (t) 


gn(t, x) dx 
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pto+h pv 0 -t+h 

+ / dt g R (t,x)dx. (39) 

J t 2 Jt — h — UQ 

ii) This case is displayed in Fig. H] arid the integral over the potential term (in the 
non-shaded area) reads 

pto+h pxo+h pto pt+h—uo 

/ / g R (t,x)dtdx — / dt g R (t,x)dx 

J t\ J X p (t) J t± J X p (t) 

rt 2 rvo —t+h rto+h pvo—t+h 

+ dt g+t, x) dx + / dt g R (t,x)dx. (40) 

J to J X p (t) J t 2 J t—h—uo 

iii) This case is displayed in Fig. 0 and the integral over the potential term (in the 
non-shaded area) reads 

rt\ rxo+h f*t\ pt+h—uo 

/ / gn{t,x) dtdx — / dt gji(t, x)dx 

Jto—hJxp(t) Jto~h J vo—t—h 

pto rt+h-uo rt 2 pVQ—t+h 

+ dt g+t,x)dx+ / dt g R (t,x)dx. (41) 

J t\ J X p (t) J to J Xp(t) 

iv) This case is displayed in Fig. Eland the integral over the potential term (in the 
non-shaded area) reads 


rt2 rxo+h 


Jx p (t) 

fto 

dt 


g R (t, x) dt dx 


rt+h—uo 


lx p (t) 


g R (t,x ) dx + 


rt2 


f^VQ—t+h 


dt 


'to 


' Xp (£) 


g R (t, x) dx. 


(42) 


For all these cases we can perform explicitly the first of the two integrals, then 
we need the explicit form of the trajectory to proceed further. We leave this for 
the direct numerical implementation since explicit expressions are straightforward but 
cumbersome. 


4. Implementation 

Here we will present the explicit implementation of the fourth order accurate algorithm 
in vacuum, the initial data and boundary conditions are also reviewed. 

4-1. Initial data 

While in the continuous Cauchy problem initial data of the second order partial 
differential equation CD are specified by ijj(t = 0,x) and dpifit = 0 ,x), a numerical 
code needs data in previous time slices; in our case two. Assuming the field evolution 
can be expanded in a Taylor series in time around t = 0 

^h,x)=+ o,x)+d t i)(o,x) (y + .. 
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where second and higher order time derivatives can be computed using the wave equation 
0 , i.e., 

3^(0, x) = 9^(0, x) - V e {x) x ) + Stm( 0, x), (44) 

we can then obtain ip(h,x) to fourth order accurate in the integration step h thus 
providing the two time levels to start the evolution algorithm. 

Note that as pointed out in Ref. [23] . this expansion requires the particle not to 
cross the line joining grid points at t — 0 and t = 2 h. One can always choose carefully 
the location of the (staggered-characteristic) grid points such that this never happens 
for particles traveling at speeds less than that of the light. 


4-2. Boundary conditions 

The boundaries of the radial coordinate x = r*, that span the space outside the black 
hole from the horizon to spacial infinity, are taken at a finite distance from the hole, x B 
and from the event horizon xa- At those boundaries we impose radiative conditions as 
purely ingoing near the horizon and purely outgoing at the outer boundary. These are 
clearly approximately true conditions which improve as we push further the boundaries. 
At the outer boundary we assume then the radiative fall off 

if(u,x) « F a (u) + ^ + ... (45) 

x or 

Evaluation of the above equation in three successive points along the constant u 
direction give us the held at the boundary as a function of the two previous levels 

ip(t + h,x B ) = ip(t, x B - h) + (46) 

x B - h) - ijj(t - h,x B - 2 h)} ^1 - + O ^ ® (^ 2 ) • 

For the inner boundary, near the event horizon xa we have instead 


which leads to 


H \ ^ , G b( v ) , G c{v) , 

Ip{v,x) « G a (v) H--- 1 -35— + ..., 


x 


X ^ 


(47) 


i>{t + h, xa) = f>(t, xa + h) + (48) 

x A + h) - ijj(t - h, x A + 2 h)] ^1 + + O + O ( h 2 ) . 

As we see, the expansions near the boundaries not only depend on the finite 
difference order, but also on the radiative condition that only applies approximately. 
This is apparent in the power dependence on the location of the numerical boundary 
and implies that pushing the boundaries further away will reduce the amplitude of the 
(spurious) reflected waves, but never completely eliminate its effects. In order to do 
that, one should use the exact treatment of the boundaries given in Ref. [T31IT3] . 
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Using Eqs. (EH)-© we can make explicit how to update the held ip using information 
on the two previous slices 


ip{t + h, x ) = —'fit — h,x) + 


i+ UI) Vt(x) 


-1 


ip{t, x — h) + f>{t, x + h) 


1 (h 


4 V 3 
4 


V £ {x — h) ip(t, x — h) + V £ (x + h) ip(t, x + h) + 16 V £ (x) fj(t, x) 


1 fh 


V £ (x - h/2) [ 1 - - f - ) V £ (x -h/2) ) x-h)+ ip(t, x)) 


1 fh 


+ V £ {x + h/2) ( ! - - ( - ) V £ (x + h/2 ) ) (ip(t, x + h) + ip(t, x)) 


■, (49) 


where (t,x) are the coordinates of the center of each cell (not a grid point), and 
g 0 = V £ (x)ip(t,x) is given explicitly in Eq. (OHlhiETTl) . 

To illustrate the numerical realization of the fourth order convergent algorithm we 
implemented the above formulae in a numerical code to compute waveforms as seen by 
an observer far away from the black holes (in this example at r* = 430M.) The initial 
distortion of the Schwarzschild black hole is produced by a Gaussian, time-symmetric 
perturbation 


ip(t = 0,r*) = A exp [— (r* — r c ) 2 /a 2 ] , (50) 

dti/>{t = 0,r*) = 0, (51) 


where A,r c and a are parameters that we have taken here as 1.0, 3.0,1.0 respectively. 

To check the convergence rate we computed waveforms at three different resolutions 
h, 2h , and Ah. In the example shown in Figs. 0and|Hl we have taken h = M/8. We 
computed the convergence rate simply as 


n = log 


ip {Ah) — ip(2h) 
ip {2 h) — ip{h) 


l l°g(2) + log|£<” l K)|/log(2), 


(52) 


where e( n )(£) represents the unknown error function of order m 1. 

The plots show the desired fourth order convergence rate on average. For 
comparison we also implemented the second order algorithm (fTTl) and displayed both 
rates together. 


5. Discussion 

Zerilli (3HJ has written down explicitly all ten held equations of General Relativity 
decomposed in tensor harmonics in the Regge-Wheeler gauge PSJ. Then derived wave 
equations for both degrees of freedom of the gravitational held, i.e. even and odd parity 
perturbations, described in terms of (gauge invariant) waveforms build up of metric 
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r*obs=430M 



Figure 7. The Waveform of a Gaussian, time-symmetric initial pulse as seen by 
an observer located at r* = 430 M. In this scale, second and fourth order evolution 
waveforms lie on top of each other. 



Figure 8. The converging power of the fourth versus the second order algorithms. 

perturbations and its first derivatives. We described here how to integrate numerically 
in the time domain these wave equation, when the perturbations are due to an orbiting 
point-like source, represented by a Dirac’s delta. In order to complete the program we 
have to be able to reconstruct the metric perturbations. This was done in Refs. PUE2, 
and applied to the radiation reaction problem in Refs. pamuu. 

While we have a good understanding of perturbations due to an orbiting particle 
around nonrotating black holes, we still need to make the equivalent progress when 
dealing with perturbations of Kerr, i.e. rotating black holes. The Teukolsky equation 
describes those (curvature) perturbations compactly as a wave equation for the 
Weyl scalars. It is straightforward to compute the energy and momentum radiated 
to infinity by the system from the scalar T 4 . Yet, we need to reconstruct the metric 
perturbations to compute, for instance the force the particle exerts on itself and that 
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accounts for the decay of orbits due to gravitational radiation (See Refs. |221 £21 for a 
review). Recently this problem has been tackled in the nonrotating limit [23 .2, • Fc> r 
the general Kerr background the problem remains open [see Ref. urn 

In any case, the starting point is to numerically solve, the Teukolsky equation with 
an orbiting particle as a source. This plays very much the role of the Zerilli and Regge- 
Wheeler equations that we just have been dealing with. A first step has been recently 
taken by generalizing the second order [TT] to fourth order [ 2 ] accurate formalism, yet 
in vacuum. The problem of generalizing the algorithms presented here for the orbiting 
particle to the Teukolsky equation remains as another open question in perturbation 
theory. 
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Appendix A. Waveforms, Source terms and initial data 


Appendix A.l. Even parity source terms 

To obtain the Zerilli equation with sources for our waveform 

[a?.-*? - vr n (r)] t) = SiZAr, t), (A.l) 


we make the following combination of Einstein equations in the Regge-Wheeler 
gauge [See Eq. (2.13) of Ref. [22] ] 

2(1 - 2 M/r) 


r( A + l)(Ar + 3 M) 


[r 2 (l - 2 M/r)d r Roo 


■(A r + M)R t 


A r 4 


ee 


(A r + 3 M) 
and End the source term is given by 


G t t + r 3 R, 


■tt 


(A. 2 ) 


(r - 2 M ) 2 (A r-r + M) A lm (r - 2 M ) 2 B £m 

(A + 1) r (A r + 3 M) * (A r + 3 M) (A + 1) 1/2 

r (A 2 r 2 + 8 rA M - r 2 A - 9 rM + 27 M 2 ) A ® 

87T —-^--—— 

(A + 1) (Ar + 3 M ) 2 

(r- 2 M) 2 V^G £m (r- 2 M)V 2 F £m 

(A + 1) (Ar + 3 M) a/A (A + 1 ) 

(r — 2 M) 3 d r A £m |C (r- 2 M)r 2 d r A ( S l 
11 (A + l)(Ar + 3M) 7F (A + l)(Ar + 3M)' 


(A.3) 
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Note that we had to reverse the sign of the (even) t 6 component of the Einstein 
equations as given by Zerilli JJ25]. Thus the source term proportional to changes 
sign. 

The source term of the Zerilli equation ra> has now the following form 


s’Zn = m fiv - m ]+gw - mi 


(A.4) 


where 


m Q U° (R-2M) [(£R) R 2 - (R-2MY 
F even (t) = — 8 tt ——————— '- Yem , 


Gevenit ) — 16 


(■ £R)nm 0 U° (R — 2 M) ±Y tm 


(XR + 3 M) (A + 1) 

(R- 2 M) U°m 0 7 rX, m (|e)|$ 

A (A + 1) 

(R-2M)U°m 0 TT ((|0 ) 2 - (sin(0 )) 2 (|$) 2 ) W £n 


A (A + 1) 


o Ytmirmo U° (R - 2 M ) 2 ((|©) 2 + (sin(0 )) 2 (|$) 2 ) 

8 (XR + 3 M) (A + 1) 

o Y lm 7T m 0 U° (R 2 A + 6 RM + 6 XR M + 3 M 2 + R 2 A 2 ) ( (±R) ) 2 

(XR + 3 M ) 2 (X + l) R 

Yimirmp U°(R- 2 M ) 2 (R 2 A + R 2 A 2 + 15 M 2 + 6 XRM) 

R 3 (A + 1) (A.R + 3 M ) 2 


(A.5) 


and where i?, 0 , $ give the trajectory of the orbiting particle in spherical coordinates 
and W £m is defined in Eq. (1 A. 28(1 . 


Appendix A.2. Even Parity Initial data 

In Refs. [EBl I24| we have introduced a conformally flat index as a combination of metric 
perturbations 


jeven _ Trim 
conf — -*-*2 








(A. 6 ) 


which is gauge invariant, and clearly vanishes for a 3-geometry that is in conformally 
flat form, with /rf m = 0, G lm = 0 and = K lm . The computation of this gauge 
invariant quantity from t/^n( r ) is most easily described in the Regge-Wheeler gauge, 
where reduces to H^ m - K lm . 













Fourth order accurate integration of black hole perturbations 


18 


Now making use of this ^ ie Regge-Wheeler gauge and the expressions 2Ti 

for Hf n and K em in terms of the waveform r fif en we can write 


jeven _ 

^ conf 


= r 


i2 „i.£m 


2 


[(\ - 3)r + 9M]M * 

i / a . 7, r\ u rTever 


(A r + 3 M)r 

[27 M 3 + 24A M 2 r + 3A(3A + l)Mr 2 + 2A 2 (A + l)r 3 ] . 


+ 


(A r + 3 M) 2 r 2 

k U°(l - 2M/r)[A(A + l)r 2 - 3M 2 ] 
(A + l)(Ar + 3M) 2 


8 [r - R] 


k U°(r — 2M) 2 

(A + l)(Ar + 3M) 


(A.7) 
d r S[r — R], 


where k = 8 Trm 0 Yt m [<d (t) , <!>(£)]. This equation is completely equivalent to the 
Hamiltonian constraint |21j . It is only written in terms of the two variables I^on? an d 
^even instead of i/| m and K (m . The conformally flat initial data that we studied in 
Refs. J21] corresponds to choose J^p = 0 and solve (with specified boundary conditions) 
for the resulting differential equation for ip^ en - 

Interestingly enough, the time derivative (represented by an overdot) of this 
equation 


'conf 


= (r - 2 M)dXen + 


a Aim , [(^ 3)r + 9 M]M - im 


(A r + 3 M)r 


-dr^ e 


+ 


[27M 3 + 24A M 2 r + 3A(3A + 1 )Mr 2 + 2A 2 (A + l)r 3 
(A r + 3 M) 2 r 2 

2k U°R[9M 2 + 2M(A - 3)r - A(A + 3)r 2 ] 

(A + l)r 2 (Ar + 3M) 2 

8 irm 0 (dY em /dt)U°(r - 2M)(A(A + l)r 2 - 3M 2 




tm 

even 


(A + l)r(Ar + 3 M) 2 

k U°r p (r - 2M) [9M 2 + 2MXr - A(A + l)r 2 


8 [r - R] 

S[r - R] 


(A + l)r(Ar + 3 M) 2 

87rmo(dYe m /dt)U 0 (r — 2 M) 2 


■d r 8[r — R] 


+ 


(A + l)(Ar + 3M) 

k U°R(r - 2M) 2 


d r 5[r — R] 


(A + l)(Ar + 3M) 


d;5[r - R] 


(A.8) 


is equivalent to the equations for the momentum constraint, i.e. when we replace the 
metric perturbations K lm , 77| m , H[ m into either the (tr) or (tip) components of the 
General Relativity constraints [2], both yield Eq. (IA.8j) above. 

By specifying conditions on J con f and / con f on an initial hypersurface E one can 
determine ffff and ip*™ which is the initial data we need to evolve Zerilli’s equation. In 
Ref. [21, we have taken 

7conf|i E 0, / con f|f E 0 . (A.9) 


to determine the “convective initial data 


ih lm f 

t even 




(A.10) 
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£m 
even I 


, d^ en {x k -,x k p [t}) 


0^#] d^en(##D 

|tE dt 


dt dx J p 




For the case of circular orbits this data takes the form 


$ 


£m | 
even I 


2m(R) \J Att/(2 £ + 1) 
A + 1 Ar + 3M 


x 


where r = l/4(y / r + y/r — 2M) 2 is the isotropic coordinate. 




£m 

even 


1 

d^ven 

L \ dt ) 

d$ 


t-E 


(A.ll) 


'yjr/r 

(A.12) 

r + y/f/r^ (R/f) e 

]f>R, 

r {y/f/r — t — lj j ( r/R,y +1 

■,r < R, 

c coordinate. 


= -imn^iZn ts > 

(A.13) 


Appendix A.3. Odd parity source term 

The Hilbert-Einstein equations in the Regge-Wheeler gauge for the odd parity sector 
(See Zerilli’sjUS] equations (C6a)-(C6c)) [Note the corrections to the source terms] 


d 2 h l 0 m d 2 h[ m 2 dh[ n 


dr 2 


drdt r dt 

_ 

(l-2M/r) v ^ATT)’ 

d 2 h[ m d 2 h £ 0 m 2 dh £ ™ 

dt 2 drdt ^ r dt 
8 ^ri(r - 2 M)Q im 

V# + 1) 

dh[ m 


+ 


AM 2(A + 1) 


ulm 

rt 0 


r — 2 M 


(A.14) 


K 


lm 


(1 - 2 M/r 


+ 2A(r - 2 M) 


dh £ n m 2 M 


(A.15) 


dr (1 — 2M/r) dt r 


h\ 


Anir 2 D em 


(A.16) 


V + 1) 

where QfJ,, Qe m and D^ m give the multipole decomposition of the energy-momentum 
tensor. 

One can use the above equations to write the metric perturbation in the Regge- 
Wheeler gauge 

CM) = ^(1 - 2M/r )d r (n/> to ) + A uffly 


(A. 17) 
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h[ m (r,t) = -- 




4mr 3 Q eri 


(A.18) 


2(1-2 M/r)^ ' A^XTT)’ 

Taking the time derivative of Eq. (IA.15D and the radial derivative of Eq. (IA.14D 
allow us to reconstruct the Regge-Wheeler equation for odd parity perturbations 


[e - a? - v^M] = s's(r,i). 


where the source is given by 


nlm _ 

^ odd 


87T (r — 2 M ) 

A a/ (A + 1) 


<9r 




d 


dt 


Im 


(A. 19) 


(A.20) 


For the Zerilli’s variable t/’ferodd the source term looks like instead 


nlm 

^ Zer,odd 


870(1 — 2 M/r) 

y/X+1 


-(1 - 2 M/r) Q lm + -^=<9,. ((1 - 2 M/r)D ern ) 


For a source term represented by a particle 
n(°) — 


(o) _ m 0 U°{t)(l - 2M/r) ang(t) 5[r - A(t)] 


Qim — 


Dtm — 


Ta/A + 1 

im 0 U°(t)(-^R) ang(f) h[r — R(t)} 
(r — 2 M)vX+T 

i m 0 R°(f) ang2(f) <5[r — R(t)} 


a/2A(A + 1) 


where 


ang (t) = 


1 fdG \ 


sin© V dt J 


) <9^ m (0, $) - sin © ( —fr ) d e Y tm {G, <f>), 


A 


rim ( 


\ dt J 


(A.21) 

(A.22) 
(A.23) 
(A.24) 

(A.25) 


ang2(t) = 


/d©\ 2 • 2 rA 

— - sm 0 

\ dt J 


x 


sin 0 


-X 




dt 

d<f> d© 


-sin0^^W^ m , 
dt dt 


(A.26) 


and -R, 0, $ define the trajectory of the orbiting particle in spherical coordinates. We 
also used Zerilli’s notation 


X* m = 2dja-cot 9)Y 


\p\VQ 


rim 


(A.27) 


w em = ( dl -cot edo - 


q 2 \ 

sin 2 9 0<f ' 


(A.28) 


The source term of the Regge-Wheeler equation (1A. 191) has now the following form 
s% = F(t) plr - R(t )] + G(i) <5[r - fi(i)], 


(A.29) 
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where 


m 


8tt m 0 U°(t) [(R - 2 M) 2 - R 2 {d t Rf 
A(A + 1 )R 


-ang(f), 


(A.30) 


G{t) = - 


87 rm 0 U°(t)R(f t R ) (J-ang(t)) 87rm 0 ang(f) 


A(A + 1) 


A(A + 1 )R 


*<* , (J t '°(*>) j t R+R2u °^ R 

+ 2 MU°(t) — RU u (t) + RU°(t) ( f n 


(A.31) 


Appendix A.f. Odd Parity Initial data 


In analogy to the ’even’ conformally flat index we define the following combination of 
metric perturbations 



(A.32) 


which is gauge invariant, and clearly vanishes for a 3-geometry that is in conformally 
flat form. With hi’f 1 = 0 in the Regge-Wheeler gauge, an initially conformally flat metric 
Icoi1\ y = 0 means 

/if* | E = 0. (A.33) 


Now taking also 




0 in the Regge-Wheeler gauge leads to 


dthT | ts = 0, 


(A.34) 


which is the Odd parity version of the CF thin-sandwich data with Uij = 0 (See 
Ref. [SHI Hn|). While / = 0 is the Wilson-Mathew j2H[ approach adapted to black 
hole initial data. Other equivalent (at our perturbative level) formulation was given 
by HUES proposing a helical Killing vector symmetry. 

From Eq. (1A.17D we see that (IA.33I) gives 


jOdd 

-*conf 



o, 


and from the definition of the waveform CD the condition h\ m = 0 gives 

f) jOdd _ n jlm \ _ rl ( L ° I 

04 con f — U — > WOdd | tE — ~^°r I —p~ J ■ 

The momentum constraint (1A. 141) then gives us the form of h 0 \ ts 


(A.35) 


(A.36) 


d 2 h 0 

oe 


+ 


2 £(I+ 1 ) 


hn 


87rm 0 U°(tjf)2M ang(fi 

(A + l) 


'*(£-&), 
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(A.37) 

where £ = r/(2M) and £ p = R(t)/(2M). 

The solutions to this equation can be written in terms of hypergeometric functions 


i tm. 


0 Its 


(A.38) 


Ci(£ p )2/i( 0, £<£ p , 

C 2 (£ p ) j/ 2 (£), £>£ P , 

[Here we have the freedom of adding y 2 for all £ as an homogeneous solution representing 
additional ’radiation content’], where 

yi( 1) = e 2 (£ - 1) _1+ ' iF 2 ([1 - £, 2 - 4 [-24 - (£ - l)- 1 ) , (A.39) 

2/2(£) = £ 2 (£ - I)- 2 "' 1 F 2 ([£ + 3 ,e + 2 ], [2 + 2 4 - (£ - l)' 1 ) . (A.40) 

In order to determine the value of the constant C(£ p ) we ensure the satisfaction 
of Eq. (I A. 3711 at r = R by equating the jump the the first derivatives of h 0 with the 
coefficient of the 5(£ — £ p ) in the source term of Eq. (IA.37I) . The result is 


Ci, 2 (£ p ) = 87rm o t/ o (i?)ang^[0(t s ),$(f s )] 


2 M y h2 


(A + 1) W 


(A.41) 


(?p) 


with W = y[ y 2 - y ' 2 yi = 21 + 1. 
For instance, for i = 2 we get 




't E 


(A.42) 


C 2 


fin 


e-i 


+ 


5(12g 3 +6g 2 +4g+3) 

3A£ 


, £>£ P - 


For £ = 3 we get 


C 1 $(2£-5/3), £ < £ p , 


ip 


1=2,m I 




c 2 


210g 3 (6g-5) 


In 


£-i 


+ 


7(360£ 4 -120£ 3 -30£ 2 -10£-3) 
2A£ 


(A.43) 
, £>£ P - 


Successive multipoles can be computed in the same way when needed. 


Appendix B. Numerical construction of up to second order derivatives 


In Sec. 13.21 we have assumed that we can Taylor expand the function g(t,x ) = 
V(x) ip(t, x) around the left and right corners of the cell (centered at (t 0 , £ 0 )) the particle 
crosses 

dg 

gn,L(t, x ) = g[t 0 , x 0 ± h\ + — [to, x 0 ± h\ (x-x 0 T h) + 

ox 

dg r, I ii /. , X , d 2 g (x — x 0 T h ) 2 

— [t 0 , xq ± h] (t — t 0 ) + [t 0 , xo ± h\ - — -h 
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Figure Bl. Case i) the particle enters and leaves to the left. The field ip is 
approximated by an expansion around the left corner of the cell using the points 
labeled by ® and around the right corner by points labeled by (g). 


d 2 g 

[to, X 0 ± h] (t - to) (x - x 0 =F h ) + 
0 2 g. (t-tA 2 




(B.l) 


So, for numerical purposes, we need to construct the derivatives appearing here in terms 
of the function evaluated in nearby grid points. The technique we use is to evaluate 
the Taylor expansion above (IB.ID in six nearby point to evaluate g and all up to second 
derivatives of g at the left/right corner of the cell the particle crosses through. 


Appendix B.l. 

Case i), the particle enters and leaves the cell on the left side. The construction of the 
derivatives is made taking into account the points labeled with the crossed circle on the 
left side, ©, and with the Xed-circle on the right, ® [See Fig. lBH ] To the required order, 
the derivatives in terms of grid points evaluations read explicitly 
For the integration on the left 


dg L 

dx 


[to,x 0 - h] = 


dg L 

dt 


[to,x 0 ~h] = 


d 2 g L 
dx 2 


[to,x 0 -h\ = 


4c?(to, x 0 - h) + g(t 0 + h, x 0 - 4 h) + g(t 0 - h, x 0 - 4 h) 

4 h 

4c/(t 0 , x 0 - 3 h) + g(t 0 - h, x 0 - 2 h) + g(t Q + h, x 0 - 2 h) 


4 h 


3c/(t 0 + h,x o — 2 h) — 3 g(t 0 — h, xq — 2 h) 

4 h 

-g(t 0 + h, x Q - 4 h) + g {t 0 - h, x 0 - 4 h) 

4 h 

g(to + h, xq - Ah) + 2 g(t 0 , x 0 - h) + g(t 0 - h, x 0 - Ah) 

4 h 2 


(B.2) 


(B.3) 


g(to - h, Xq - 2 h) + </(t 0 + h, Xq - 2 h) + 2c/(to, x 0 - 3 h) 

4h 2 
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d 2 g L r , ,, g(t 0 -h,x Q - 4 h) - g(t 0 -h,x 0 - 2 h) 

-[t 0 , x 0 - h\ = - 


dtdx 


Ah 2 

-ff(to + h, x 0 - Ah) + g(t 0 + h, x 0 - 2 h) 
Ah 2 


d 2 gL u ,, 3g(t 0 - h, x 0 - 2 h) - 2 g(t 0 , x 0 - h) 

-Rn^o - © = - 


c)t 2 


(B.4) 


(B.5) 


4/z 2 

g(t 0 + h, x 0 - Ah) + g(t 0 + h, x 0 - 2h) 

Ah 2 

g(t 0 - h,x 0 - Ah) + 3^(t 0 + h, x 0 - 2 h) - 6 g(t 0 , x 0 - 3 h) 


Ah 2 


and for the integration on the right 


(B.6) 


9g H 

dx 


dgn 

dt 


[t 0 ,x 0 + h} = 


[to, xq + h] = 


d 2 g R 
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d 2 g R 

dtdx 


d 2 g R 
dt 2 


[t 0 ,x 0 + h\ = 


-gjfo + h,x o + 2 h) + gjto - h, x 0 ) 

4/i 

<7(t 0 + h, x 0 ) - (7 (t 0 - h, x 0 + 2h) 

+ Ah 

g(t 0 + h, x 0 ) - g(t 0 -h,x 0 + 2h) 

Ah 

-g (to - h, x 0 ) + g(to + h, x 0 + 2h) 

Ah 

-g(tp + h, xp + 2 h) - 2g(t 0 , x 0 + h) + g(t 0 - h, x 0 ) 
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g(tp + h, xq) — g(tp — h,x o + 2 h) + 2 g(tp, Xp + Sh) 
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(B.8) 


Ah 2 


(B.9) 


[t 0 , x 0 + h] = 


[to, Xq + h] — 


(t 0 ~ h, Xp) + (tp + h, Xp) 

Ah 2 

-g(t 0 + h, xp + 2 h) + p(t 0 - h, x 0 + 2h) 

4/r 2 

3ff (t 0 - h, x 0 + 2h) - 6g(t 0 , x 0 + h) - 2g(t 0 , x 0 + 3h) 

Ah 2 

g(tp — h, Xq) + g(tp + h, Xp) + 3 <7 (to + h, Xp + 2h) 


(B.10) 


4h 2 


(B.ll) 


Appendix B.2. 

Case ii), the particle enters the cell on the right side and leaves it on the left side. 
The construction of the derivatives is made taking into account the points labeled with 
the crossed circle on the left side, ©, and with the Xed-circle on the right, <g) [See 
Fig. IB21 ] Note that if we had chosen the alternative point g(t 0 + h, Xp — Ah) instead 
of g(tp — h,x o — Ah) to the left, the set of equations do not produce a solution. The 
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x0-4h xO-3h xO-2h xO-h xO xO+h x0+2h xO+3h x0+4h 

Figure B2. Case ii) the particle enters the cell on the right and leaves it on the left. 
The field ijj is approximated by an expansion around the left corner of the cell using 
the points labeled by ® and around the right corner by points labeled by (g>. 


same applies to the right of the particle if we had taken g(t 0 — h,x 0 + 4 h) instead of 
g(t 0 + h, x 0 + 4 h) as the sixth grid point. 

To the required order, the derivatives in terms of grid points evaluations read 
explicitly for the integration on the left 


dg L , ,, — 2 g(t 0 , x 0 - 3 h) + 2 g(t 0 , x 0 - h) - 2 g(t 0 - h, x 0 - 2 h) 

— [t 0 ,x 0 -h] = 


dg L 

dt 


[t 0 ,x 0 - h\ = 


d 2 g L 


[t 0 ,x 0 - h] = 




Ah 





- h,x 0 ) +g(t 0 - 

- h, x 0 

— Ah) 



1 

Ah 


1 



-2 g(t 0 

— h, xo — 2 h) — 

O 

1 

h, x 0 ) + 2 < 7 (t 0 , 

xo ~ h) 




Ah 




-2 g(t 0 , 

Xq - 3 h) + <7 (t 0 

- h,x o - Ah) + 2g(t c 

i + h, Xo 

— 2 h) 



Ah 



-2 g(t 0 

— h, xo — 2 h) + 

1 

O 

-40 

h,x 0 ) + g (t 0 — 

h, x 0 - 

Ah) 


4 h 2 


(B.12) 

(B.13) 

(B.14) 


dx 2 

d 2 g L r n _ -g(to - h, x 0 ) + 2 g(t 0 , x 0 - h) - 2 g(t 0 , x 0 - 3 h) + g(t 0 -h,x 0 - 4 h) 

[to, xq — n\ — 


dtdx 

SV 

dt 2 


A1F 


[t 0 ,x 0 - h] = 


gjtp -h,x 0 - Ah) - Ag(t 0 , x 0 - h) + g(t 0 - h, x 0 ) 

Ah 2 

-Ag(t 0 , x 0 - 3 h) + Ag(t 0 + h,x 0 - 2h) + 2 g(t 0 -h,x 0 - 2h) 

Ah 2 


(B.15) 


, (B.16) 


and for the integration on the right 


dg R 

dx 


[to, Xq + h] — 


dg R 

dt 


[to, Xq + h] — 


- 2 p(t 0 , Xq + 3 h) + 2 g(t 0 , x 0 + h) - 2 g(t 0 + h,x 0 + 2 h) 

Ah 

g(t 0 + h, x 0 ) + g(t 0 + h, x 0 + Ah) 

+ Ah 

9 (to + h, x 0 ) - 2<7 (t 0 , x 0 + h) - 2 g(t 0 - h, x 0 + 2 h) 

Ah 


(B.17) 
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Figure B3. Case iii) the particle enters the cell on the left and leaves it on the right. 
The field i/j is approximated by an expansion around the left corner of the cell using 
the points labeled by © and around the right corner by points labeled by ® . 


d 2 g R r 

-JL[t 0 ,x 0 + h] 
d 2 g R 

4l?[to,io + h] 


2 g(to + h, x 0 + 2 h) — g(t 0 + h, x 0 + Ah) + 2 g(t 0 , x 0 + 3 h) 

I--4ft-• (B. 18) 

-2g(t 0 + h, Xq + 2 h) + g(t 0 + h, x 0 ) + g(t 0 + h, x 0 + 4 h) 

- Th? -' (B19) 

-git o + h,x 0 + 4 h) - 2g(t 0 , x 0 + h) + 2 g(t 0 , x 0 + 3 h) + g(t 0 + h, x 0 ) 


4 h 2 

4gjt 0 - h, x 0 + 2/t) - Agjt 0 , x 0 + h) + 2g(t 0 + h, x 0 + 2/t) 

4/t 2 

g(to + h, Xq) + ^(t 0 + h, Xp + 4/t) — 4(yf(t 0 , % + 3/t) 

+ 4^2 • 


5 

(B.20) 


(B.21) 


Appendix B.3. 

Case iii), the particle enters the cell on the left side and leaves it on the right side. 
The construction of the derivatives is made taking into account the points labeled with 
the crossed circle on the left side, ©, and with the Xed-circle on the right, ® [See 
Fig. IB31 ] Note that if we had chosen the alternative point g(tp — h,x o — 4 h) instead 
of g(t 0 + h,xo — 4 h) to the left, the set of equations do not produce a solution. The 
same applies to the right of the particle if we had taken g(t 0 + h,x 0 + Ah) instead of 
g(tp — h,x o + Ah) as the sixth grid point. 

To the required order, the derivatives in terms of grid points evaluations read 
explicitly for the integration on the left 


dg L u _ -2 g(t 0 , x 0 - 3 h) + 2g(t 0 , x 0 - h) - 2g(t 0 + h,x 0 - 2h) 

g(t 0 + h, x 0 ) + g(t 0 + h, x 0 - Ah) 

+ Ah 

dg L r ,, g(to + h, x Q ) - 2g(t 0 , x 0 - h) - 2g(t 0 -h,x 0 - 2h) 

— [t 0 ,x 0 - h\ - 


(B.22) 


Ah 
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2 g(t 0 + h, x 0 - 2 h) - g(t 0 + h,x 0 - Ah) + 2g(t 0 , x 0 - 3h) 
+ 4 h 

d 2 g L , . 1 _ —2g(t 0 + h, x 0 - 2/i) + g(f 0 + h, x 0 ) + g(t 0 + h, x 0 - Ah) 

'FOj^O — ,l \ — - 


dx 2 


Ah 2 


, (B.23) 

(B.24) 


d 2 g L r ,, ~g(t 0 + h,xo- Ah) - 2g(t 0 , x 0 - h) + 2p(t 0 , x 0 - 3 h) + p(t 0 + h, x 0 ) 

ho, x 0 - h\ = - 


dtdx 

d 2 g L 

dt 2 


[t 0 ,x 0 - h\ = 


Ah 2 


4< 7 (t 0 - h, x 0 - 2h) - Ag(t 0 , x 0 - h) + 2g(t 0 + h, x 0 - 2h) 

Ah 2 

< 7 (t 0 + h, x 0 ) + < 7 (t 0 + h, x 0 - Ah) - Ag(t 0 , x 0 - 3h) 


Ah 2 


(B.25) 


(B.26) 


and for the integration on the right 


dg R -2g(t 0 ,x 0 + 3h) + 2g(t 0 ,x 0 + h) - 2g(t 0 - h,x 0 + 2h) 

-[t 0 ,x 0 + n\ — - 


dx 


Ah 

g(t 0 - h, x 0 ) + g(t 0 -h,x 0 + Ah) 
Ah 


(B.27) 


dg R -2 g(t 0 - h, x 0 + 2 h) - g(t 0 - h, x 0 ) + 2g(t 0 , x 0 + h) - 2g(t 0 , x 0 + 3 h) 

■[to, xq + h\ = - 


dt 


Ah 

<7(to - h, xq + Ah) + 2 g(t 0 + h, x 0 + 2 h) 
Ah : 


d 2 g Ru —2g(t 0 — h,x 0 + 2h) + g(t 0 — h, x 0 ) + g(t 0 — h,x 0 + Ah) 


' [t 0 , Xq + h\ = 


Ah 2 


(B.28) 

(B.29) 


dx 2 

d 2 g R -g(t 0 -h,x 0 ) + 2g(to,x 0 + h)-2g(to,x 0 + 3h)+g(to-h,x 0 + Ah) 

ho, Xq + n\ — 


dtdx 

9 2 g« 

dt 2 


Ah 2 


(B.30) 


[to, Xq + h\ — 


g(t 0 - h,x 0 + Ah) - Ag(t 0 , x 0 + h) + g(t 0 - h, x 0 ) - Ag(t 0 , x 0 + 3 h) 


Ah 2 

4< 7 (t 0 + h, Xq + 2 h) + 2g(t 0 - h, x 0 + 2 h) 
Ah 2 


(B.31) 


Appendix B.f. 

Case iv), the particle enters and leaves the cell on the right side. The construction of 
the derivatives is made taking into account the points labeled with the crossed circle on 
the left side, ©, and with the Xed-circle on the right, <g) [See Fig. IB4I ] To the required 
order, the derivatives in terms of grid points evaluations read explicitly 
For the integration on the left 


dg Lf , n _ -g(t 0 + h,x 0 - 2h) + g(t 0 - h,x 0 ) + g(t 0 + h,x 0 ) - g(t 0 - h,x 0 - 2h) 

-[to,XQ — n\— - 


dx 


Ah 


(B.32) 


dg L r . n _ g(to + h,x 0 ) - g(t 0 - h,x 0 - 2h) - g(t 0 - h,x 0 ) + g(t 0 + h,x 0 - 2h) 
- W [ t o^o-h\- — : 


(B.33) 
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x0-4h xO-3h x0-2h xO-h 


xO xO+h x0+2h xO+3h x0+4h 


Figure B4. Case iv) the particle enters and leaves the cell on the right. The field 
•0 is approximated by an expansion around the left corner of the cell using the points 
labeled by © and around the right corner by points labeled by 


d 2 g L u n _ -g(t 0 + h, x 0 - 2h) - 2 g(t 0 , x 0 - h) + g(t 0 - h, x 0 ) 

dx 2 ^ 0, X ° J “ Ah 2 

g(t 0 +h, x 0 ) - g(t 0 - h,x 0 -21i)+ 2g(t 0 ,x 0 - 3h) 

+ -^-. < B - 34 > 

d 2 g L u n _ —g(t 0 ~ h, x 0 ) + g(t 0 + h, x 0 ) - g(t Q + h, x 0 - 2h) + g(t 0 - h, x 0 - 2h) 
© 0 ) x o ~ ll \ — - 


dtdx 

dpgi 

dt 2 


Ah 2 


[to,x o - h\ = 


3 g(t 0 -h,x o - 2 h) - 6 g(t 0 , x 0 - h) - 2g(t 0 , x 0 - 3 h) 

Ah 2 

g(to — h, Xq) + g(to + h, Xq) + 3 g(t 0 + h, x$ — 2h) 


Ah 2 


and for the integration on the right 


(B.35) 


(B.36) 


dg R 

dx 


dg R 

dt 


[t 0 ,x 0 + h\ = 


Ag{t 0 , x 0 + h) + g(t 0 + h, x 0 + Ah) + g(t 0 ~h,x 0 + Ah) 

Ah 

Ag{t 0 , x 0 + 3 h) + ( 7 (t 0 - h, x 0 + 2 h) + g(t 0 + h, x 0 + 2h) 


Ah 


(B.37) 


[t 0 ,x 0 + h} = 


d 2 g R 
dx 2 


d 2 g R 

dtdx 


[to, xo + h\ = 


3g(t 0 + h, xo + 2h) — 3g(t 0 — h,x o + 2 h) 

Ah 

-g (to + h,x o + Ah) + g {to- h, x 0 + Ah) 

Ah 

g(t 0 + h,x o + Ah) + 2g(t 0 , Xq + h) + g{to — h, Xq + Ah) 

Ah 2 

g{to — h,x o + 2 h) + g(to + h, Xq + 2 h) + 2g{t 0 , Xq + 3 h) 


(B.38) 


[t 0 , Xo + h\ = 


Ah 2 


g{t 0 -h,x o + Ah) - g{t 0 -h,x 0 + 2h) 
Ah 2 

—g(to + h, x o + Ah) + g(to + h, Xq + 2 h) 
Ah 2 : 


(B.39) 


(B.40) 
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d 2 g R 

dt 2 


[to, xq + h] 


3g(t 0 —h,x 0 + 2 h) - 2g(t 0 , x 0 + h) 


Ah 2 

~\~g{to + h, Xo + Ah) + g{to + h, xq + 2 h) 

Ah 2 

A~g(t o — h, Xq + Ah) + 3g(fo + h, Xq + 2h) — 6g(fo, %o + 3 h) 


Ah 2 


(B.41) 
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